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ABSTRACT 


A method for numerically solving the three-dimensional unsteady 
Euler equations using flux vt&edi splitting is developed. The equations 
are cast in curvilinear coordinates and a finite volume discretisation 
is used. An explicit upwind second-order predictor-corrector scheme 
is used to solve the discretized equations. The scheme is stable for a 
CFL number of 2 and local time stepping is used to accelerate convergence 
for steady-state problems. Characteristic variable boundary conditions 
are developed and used in the far field and at surfaces. No additional 
dissipation terms are included in the scheme. Numerical results are 
compared with results from an existing three-dimensional Euler code 
and experimental data. 



I. INTRODUCTION 


In support of NASA's interest in the use of prop-fans as a propul- 
sion device, a computational method was recently developed for numerically 
solving the flowfield about a swept, tapered, supercritical tiring with a 
propeller producing thrust and swirl (Ref. 1). The equations solved 
were essentially Euler equations with body force terms included to 
simulate the propeller. The Euler solver used was an extension of the 
method of Jameson, et. al. (Refs. 2 and 3) referred toasFLQ57. Although 
good agreement was obtained with experimental data (Ref. 1), some 
difficulty in convergence, particularly the energy equation, was encountered 
using the FL057 central-difference scheme. Recently, convergence 
difficulties were also reported by Swafford (Ref. 4) in solving a set 
of hyperbolic equations with source terms (which can be considered 
similar to the force terms included in the Euler equations in Ref. 1) 
using a similar central-difference scheme. An upwind scheme eliminated 
the convergence problems in Ref. 4. Moreover, any additional smoothing 
added to the upwind scheme was found to be detrimental with regard to 
convergence. (Additional smoothing was found to always be necessary 
in using the central-difference scheme in Ref. 4,) Therefore, it 
seemed appropriate to consider an upwind scheme for solving the three- 
dimensional Euler equations. 

The flux-vector-split form of the equations used is developed in 
the following section. The motivation and background of using splitting 
is given in the literature (see, for example, Steger and Warming (Ref. 5)) 
and is not repeated here. Moreover, many of the matrices needed are 
given in the literature (see Refs. 6 and 7) ; however, all equations 



and matrices needed are developed and included in Section II for clarity 
and completeness. Formulation of the equations for numerical solution 
and the algorithm used to solve the discretized equations are discussed 
in Section 111. 'tsh&raet*SKistic variable boundary conditions used 

in the far field and at surfaces are developed in Section IV. Numerical 
results, including comparisons with FL057 solutions and experimental 
data, are presented and discussed in Section V. 
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, II. SPLITTING 

The conservation law vector form of the -Easier equations in 
Cartesian coordinates x,y, and z are 


3h 


JL§. 4 . _ p- 

3t ^ 3x + 3y + 3z " G 


( 2 . 1 ) 


where 


q ** Ip» Pu, PV, pw, e] 

2 T 

f = [pu, pu + p, puv, puw, u(e -*• p)] 

2 T 

g = [pv, puv, pv + p, pvw, v(e 4= p) ] 

2 T 

h “ [pw, puw, pvw, pw + p, w(e -iir p) J 
p = (y " 1) [« - %p(u 2 -f V 2 + w 2 } J 


Using curvilinear coordinates defined as 

£ = £(x, y,z) 
n " n(x,y ,z) 

C “ S(x,y,z) 

T ■= t 


it is straightforward to transform Eq. (2.1) to 

iS + !£ ' + ± G - + ' iH . 0 

St 35 an 3 ? 


( 2 . 2 ) 


Q « J[p, pu, pv, pw, e] 


where 



F « J[pU, puU + £ x p, pvll + CyP» PwU + £; z P, U(e + p) ] 

T 

G *= J[pV, puV + n P, pvV + n p, pwV -I- tip, V(e -f p) ] 

, a y « 

H = J[pW, puW + C X P, pvW + E y P, pwW + C z p, W(e + p) ] T 

J " V Vs ~ Vs* " VVs " z nV + VVs ‘ y nV 


- J (y n * c - « n y c ) 


J ~ 1( Vs - VV 


5 x " J ” <V n ' Vn } 
S " J ~ 1( Vc - Vs } 


ORIGINAL FACE KZ 
OF POOR QUAL.!T¥ 


n y " J ( V'C ' Vt> 
S = J " 1(z e x n - Vn> 
5 Z = J ~ 1( Vr, ' y nV 


J “ ( Vs " Vs } 


J (x..y - y r x ) 


u " V + * y v + V 


V a i)u + nv+nv 
x y z 


W » ( u + ( v + 5 w 
x y J z 


The strong conservation .law of the Euler equations in curvilinear 
coordinates (Eq. (2.2)) can be written in the quasilinear form 


|| + A || + B || + C || - 


0 


(2.3) 



V? 


where the matrices A, B, and C are given by 

8F 
SQ 

3G 
3Q 

3H 
9Q 


A 

B 

C 
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Carrying out these operations, one obtains a matrix K 


K 


kj - ue R 


k <j> - V0. 
y k 


v - we k 


k x ( 2 -- y)u + e R 


k v - k (y - l)u 
x y 


k w ~ k (y - l)u 
x z 


( 24 * - xe -)9, k ( x ® - - (y -■ Due. 

P iv P 


y 


k u - k (y - l)v 
y X 


k (2 ~ y)v + 0. 
y k 


k w ~ k (y - l)v 
y z 

k y C x - - 0 - Cy -Dve k 


k 2 u - k x (y - l)w 


k v - k (y - l)w 
z y 


k,(2 - y)w + 0 k 


k (3-°- - <t > ) - (y - l)w0 

Z p r 


k x (y - 1) 


k y (y - 1) 


k 2 (y - 1) 


ve t 


(2.4) 


where 


4 > ■- (u 2 + + w 2 > 

0, « k u + k v + k w 
■k x y a 


and matrices A, B, and C tire given by the matrix k depending on whether 


r‘sl 
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k in Eq. (2.4) is n s or £. That is 

ic = A for k ** £ 

K = B for k = n (2.5) 

K *= C for k ° c 

The eigensystem of the matrices A, B, and C is important to 
achieve -he intented splitting. However, there are few zero elements 
in Eq. (z.4) and the eigenvalues of A, B, and C are difficult to 
determine using Eq. (2.4) directly. Hence, consider the nonconservative 
vector form of the Euler equations in curvilinear coordinates 


la + a la. + b la + c la. _ 0 

9t 8£ 3n 8r, 


( 2 . 6 ) 


where 

q *=■ J[p, u, v, w, p] 
Note that Eq. (2.3) can be written as 


T 


M la am + BM |~ 3 + CM -I 3 - = & 
8t 8£ 8'0 8C 


(2.7) 


where M is the matrix Multiply Eq. (2.7) on the left by M ^ 

«q 

to obtain 

I I 3 -I- M -1 AM I 3 + M _1 BM I 3 + M -1 CM | 3 - 0 
3 T 3£ 3n 3C 

where I is the identity matrix. Then from Eqs. (2*6) and (2.7) 

a - M" 1 * AM 


b « M BM 
c •= M~ 1 CM 


( 2 . 8 ) 
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Therefore, the matrix a is similar to A, b is similar to B, and c is 
similar to C. Because similar matrices have the same eigenvalues 
(Ref. 8), the eigenvalues of A, B, and C are Known by determining the 
eigenvalues of a, - and e'=. The matrices a, b, and c are more simple 
to handle than matrices A, B, and C because they contain several zero 
elements as shown below. 

The matrices a, b, and c are now determined. The matrix M gi’ en 


by 


IS 

8q 



Y-l 


PU 


pv 


0 0 

0 0 

0 0 

P o 



( 2 . 9 ) 


The inverse of this matrix is 



<J> . -u(y~1) -v(y-1) 


0 

0 


0 


( 2 . 10 ) 



-w(y~1) 


(y-D 
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Then by Eqs. (2.8) the matrices a, b, and c are given by 


r 


K « 


9 k pk x 


P k 


? 2 
0 k pc k pc 
* y 


Pk 


k z pc 


x 

P 

k 

.J l 

P 


k 


where 


k ~ a for k - £ 

k ~ b for k = n 

k = c for k ~ £ 


The c in the pc terms in Eq. (2.11) is the speed of sound 
to be confused with the matrix e. The $ appearing in Eqs. 
and (2.10) is again 

4> = (u 2 + + w 2 ) 

just as in Eqs. (2.4). 

From Eq. (2.11) the eigenvalues of matrices a, b» and 
easily determined. The eigenvalues are 

x t‘'4" x k*V + S’ + Y’ °k 


x t ” \ + c l vk -l 


K " \ - 


( 2 . 11 ) 


( 2 . 12 ) 

and is not 
(2,9) 


c are 


(2.13) 


/ 



where A*, A^, A^, A^, and A^ are the eigenvalues of a for k = £» b for 
k = n* and c for k - £. Correspondingly, the eigenvectors are 


where 


x, « [k , 0, k , -k , 0] T 

1 x* * z’ y 

x, *» [lc , -k , 0, k , 0] T 

2 y z’ ’ x’ 

x~ ■= [lc , k , -k , 0, P 1 '*' 

3 z* y* x’ ’ 

i r p r ~ .t 

x, = l— , K , k , k , pc] 

4 J2 C x y z 

x, = — [— , -k , -k , -k , pc]* 

5 c x y z 


k = 


x 


x Vk 


k y ° Tvkl 


2 J t 2 ^ 1 2 ‘ ^ 

lc + k + k r 

x y 2 


k « 


TW 


( 2 . 14 ) 


( 2 . 15 ) 


Sufficient equations have now been developed to carry out the 
splitting. The integral conservation law form of the Euler equations 
will be formulated and discretized for numerical solution in the 
following section. This formulation requires evaluation of the 
vectors F, G, and H at faces of the finite volumes. By splitting 
the vectors F, G, and H into the sum of separate vectors, each having 
an eigenvalue as a coefficient, the evaluation of each separate vector 
by extrapolating from the appropriate direction indicated by the sign 
of its eigenvalue can be carried out. 
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The vectors F, G, and H are homogeneous functions degree one 
in Q; therefore, by Euler's Theorem (Ref. 9) 


K = KQ 


•’2.16) 


where K = F and K = A for k - £, K = G and K = B for k = n, and 
K = H and K ■= C for k = £. The matrices K can be written 


K » T.A, T. 
k. k k 


(2.17) 


where is the diagonal matrix whose diagonal element? are the eigen- 
values of k (cr, also, K) given by Eq. (2.13) The matrices k given 
by Eq. (2.12) can be written 


k >» p. a, pr 1 

k k k 


(2.18) 


where the columns of P^ are the exgem ector s of k corresponding to the 
respective eigenvalues. From Eqs. (2.8) 


K ■» MkM 


(2.19) 


Using Eq. ^2.18) in Eq. (2.19) 


k = mp. a, prV 1 

k k k 


( 2 . 20 ) 


Then from Eqs. (2.17) and (2.20) 


T, -- MP. 
k k 


(2.21a) 


-1 -1 -1 

T, = P. M 

k k 


(2.21b) 


The matrices M and M ' are given by Eqs. (2.9) and (2.10). 


matrix P. is 
lc 
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where 



Rather than invert Pj directly, the matrix can be determined easily 
from the left eigenvectors of k. Vhe left eigenvectors of k are the rows 






1st' 



where 


OF POOR QUALITY 



& 



;• 

/2 pi ' ; 



-1 

The matrices and 

can now be determined using Eqs, (2.21). 

V 

These 

matrices are 



j. 


k 

k 

k 



X 

y 

z 






■ » 


uk 

X 

uk - pk 

y 2 

uk + pk 
z y 






J 

T ** 
k 

vk + pk 
x z 

vk 

y 

vk z - pk x 

1' 


wk - pk 
x y 

wk + pk 

y 11 

Wk^ 



-JiL. ^ .j. p(vk - wk ) 
y-.l x K z y 

— k 4- p (wk — 
Y-l y x 

uk ) • '^r k 4* p(uk - vk ) 

z * y 

- « 
* V? 

a 

a 

. 


a(u 4* ck ) 

A 

a(u - ck ) 

X 






j 

a(v + cky) 

Ct(v - ck y ) 

(2.24a) 

.* 

a(w + ck,,) 

a(w - ck ) 

p. 


* ;; 

, l( ±g! + c ; k ) 

- cO.) 
y~l k 




2 > 



OWGINAV- 
OF 


-1 


k x (l- ~2*) + ~(wk -vk z ) 


k ( 1- 4) ( 


k z (l- -|) + “(vk x -uk ) 


3<* - c6 k ) 


+ cO k ) 


-k i + k ^ 

yp X c 2 

U + k 

x p y c 2 

2 2 
c 


e[ck z - w( y-l) ] 


-f3[ck z + w(Y~l)l 


j; 

x 2 
c 


v(.Yrii 


^ +k x 2 

c 


-ki + k 


ZP y c 2 


k- + k — 


yp 


y 2 

J c 


-k -- + k„ .. 

Kp s 2 

K c 


V(1-J0 


3fck - u(y- 1)1 Pick - v(y - 1) 1 

X y 


-3[ck x + u(v-l)I -P [cky + v(y~1)1 


-k^=| 
x 2 
c 

_kl=l 
y 2 
* c 

_k I=i 

k z 2 
c 

P(Y-l) 


e(Y-D 


(2.24b) 


where. 


Y-l , 2 , 2 2. 

-Jy- (ll + V + W ) 


0, * k u + k v + lc w 

k x y z 


_JL_ 
& c 


1 

Jl pc 


Using Eqs. (2.16) and (2.17) 



The diagonal matrix 
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(2.26) 


can be written 


1 4 % 

A. I. _ „ + A, I, + A,‘ I r 
k 1,2,3 k 4 &. 5i 


(2.27) 


where 1^ ^ 3 is a matrix which has unity as the: first three diagonal 
elements and all other elements zero, has unity as the fourth diagonal 
element and all other elements zero, and has unity as the fifth 
diagonal element and all other elements zero', fete that only three 
terms are needed on the right hand side of Eq, (2,27) instead of five 
because the first three eigenvalues are the same as shown by Eqs. (2.13). 
Using Eqs. (2.25) and (2.27), the split form, of K is 

K ■ >t T U h.2.3 <*« + \ h \ h ^0 < 2 - 2S > 

K = K x + K 2 + K 3 


or 



where 


ORIGIN At. V-Aciii !j 

OF POOR QUALITY 


mt* i \ h.2,3 T k 1 « 

k 2 ■ X J \ h 1 

K 3 " l k T l« h T k l 1 


Everything is now available to carry out these operations and determine 
the split form of the vector K. The operations described by Eq. 

(2.28) give 

K ■*» K x + K 2 + K 3 (2.29) 


where 


K = 

F 

for 

k =■ 

K. = 

G 

for 

k = 

K = 

H 

for 

k = 


K * •?— 

1 k y 


pu 


pv 


pw 

2 2 2 
kP (u + V + w ) 



\ 
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pu + pck 


K 2 - X* -i PV + pck, 


pw -f pck 


e + p + pc0. 


P« - pck 


5 T 

K_ » A. -r— pv - pck 
3 k 2Y * y 


pw - pck 


e + p - pc0. 

• If 


X k = A k = A k" U X U + V + k 2 W " °k 


4 - 0 k + c | Vk| 

A k“ \~ c l Vk l 

|Vk| - (k* + k* + k*) % 


1 | Vk| 


e, » k u + k v + k w 
k x y z 


This is the same form of the equations as presented by Reklis and 


Thomas (Ref. 10). 
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The discretized integral fox® of Eq. (2.2) in computational 
space for a cell with center denoted as i,j,k is (see, for example, 
Ref. 11) 


A3 

At 


ASAqAC + (F i+% 


- * ± .*> AnAt + < G j+ 


j- 

'z 



or 


+ (H 


k+’a 


H, , )A£Aq = 0 

K-'g 


0 . 1 ) 


AS 

At 



+ — 1 + 
An 


6 H, 

k 

A l 


= 0 


(3.2) 


The central difference operator notation to Eq. (3.2) indicates the 
flux vectors are evaluated at cell faces to this finite volume formu- 
lation as illustrated in Fig. 1. 

The numerical scheme used is a finite volume version of the 
second-order upwind scheme of Warming and Beam (Ref. 12). The 
present scheme is an extension of that used by Deese (Ref. 13) for two- 
dimensional flow. To illustrate this scheme consider the simple model 
equation 


It MrA „ 0 

3t 3x 


(3.3) 


where. F(u) *= au; and a is a constant taken as greater than zero 

for illustration. The second-order upwind scheme of Warming and Beam 


is 


-n+1 




At 


Ax 


(3.4a) 
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V^F n yp n+ ^ 

n+1 , , n . -n-ML At i At . i 

u ± " ^ (u i + u i 3 " t ~ax t — sr 


(3.4b) 


where 


S^i*.a^ n .ft5'p n -» F(u n ) - F(u n ,) 

■ i i-i u r ^ i-r 


-n+1 , , -n+1. 

F i “ F(u i > 


The finite volume form of Eq. (3.3) corresponding to Eq„ (3.2) ±3 


A 6F 4 

Au , i „ 

~T~ + -7 “ 0 

At Ax 


(3.5) 


By using the one- point upwind extrapolations, u for u. , and u , 

1 1+7 1-Jl 

for u , , the predictor step for Eq. (3.5) is 
X-'t 

VF n 

-n+1 n aAt, n n . n .. 1 , ~ 

u i c u i - ir (u i “ Vi* “ u i - (3 * 6) 


which is the same as Eq. (3.4a). By using the two-point upwind 

extrapolations 2u. - u , for u., f and 2u. . - u. 0 for u. , , the 
i i-1 i+% i-1 i-2 

finite volume corrector step is 


n+1 n aAt r .„ n n . n n . , ,-n+l -n+1*, 

u i = u i - m [(2u i “ u i-i> ~ C2u i-i ” u i-2 } + (u i - u i-i )3 


n+1 n aAt r/ n „n , n , , , n n At i 

U i = U i “ 2Ax U i - 2u i-l + U i-2> + (u i “ u i-l )] - ~2 Ax 


Using Eq. (3.4a), Eq. (3.7) can be written 


n+1 , / n . -nih At i At VF 

u . = k\u. + u. ) — - -= — 

i i 1 l Ax 2 Ax 


(3.8) 
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which is the same form as Warming and Beam's corrector step given by 
Eq. (3.4b). The stability, dissipation, and dispersion analyses of 
Ref. 12, tawfdfere, are applicable to the present finite volume scheme. 
In addition, it is not difficult to show that this scheme is also 
consistent. 

The scheme for the complete three-dimensional unsteady equation 
is: 

Predictor: 


+ ^.J-rt.k’ - V'ol.i-H, k> + 'U<«i.J,krt> 

- »< ik-^i »- 9 > 


where the subscript t corresponds to one part of the split vector in 
Eq. (2.29); and 

**11 n 14 3 

Q., t , , . = Q. . , if the corresponding X , X or X eigenvalue 
i+'ziJ » *t r,j,K i, t, £, 

evaluated at i+%,j,fc is > 0 

Q^., < i " Q?. , . . if the corresponding X^.xjl, or X^ eigenvalue 

evaluated at i-Ra,j,k is < 0 

and similarly for Q P , ... The same is done for (f 1 , . , . and 

Q 1 ? i 7 t except x\ X , and X eigenvalues are interrogated; and 

r > j -•£ > nun 

*“n *”n 14 5 

similarly, for Q , . and Q , , except X X , and X 

eigenvalues are interrogated. 





Corrector; 
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"ilk ’ <?j,k - -V* **d V3«,.j.k> - 


An ah a« a« 

+ V Q i,j+^,l^ + VViW ~ V Q i,j,kV 

+ vCL,^ - vCL, k > + G *Kku.> - v?*U, k > 




(3.10) 


~n*t\L 

where the Q , . . , etc, variables are determined in the same way 

j *k 

as the 5 “ . . , etc, variables as described below Eq. (3.9) except 

— n A n 

Q , etc, are used in place of Q . , etc. The Q . , , etc, are 

3* > i * * * J > > J 9 ^ 


determined by 


.1 ,4 


Q i+W.k ■ 2Q i.j.k ‘ Vl.j.k i£ tl,e A e • V ° r \ 

eigenvalue evaluated at i+^,j,k 


*n „.n ^n 

Q i+%,j,k " 2Q i+l “ Q i+2,j,k 


is > 0 


1 A 5 

if the corresponding X A or X 

^ s s» 

eigenvalue evaluated at i+^,j,k 


is < 0 


mi 

and similarly for Q. , , , . Again, the same is dorse for Q. , , . 

and Q 1 ? . . . by interrogating the appropriate r| or 5. eigenvalues. 
i»J 

Although the algorithm given by Eqs. (3.9) and (3.10) is an 
extension to three dimensions of that used in Ref. 13 for two dimensions, 
the interrogation of eigenvalues differs significantly. Three methods 
were tried in Ref. 13 to handle the computation of flow variables at 
cell faces when eigenvalues were of different sign on either side of 
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the cell face. The approach taken in the present scheme is to compute 
eigenvalues at cell faces rather than at cell centers as illustrated in 
Fig. 2 on a £ « constant ^lai^c The difficulty associated with eigenvalues 
changing sign is thus eliminated. This seems a natural approach because 
the use of Eq. (2.29) in the finite volume discretized equations (Eq. (3.2)), 
requires that eigenvalues be known at cell faces. The eigenvalues are 
computed by averaging the information on either side of a cell face that 
is necessary for their computation according to Eq. (2.13). Then, 
depending on the sign of the eigenvalues, the information necessary to 
compute the remaining terms in the split vectors given by Eq. (2.29) for 
use in the algorithm given by Eqs. (3.9) and (3.10), is determined by 
extrapolation from the appropriate direction. 

The time step At., t is determined from 

(3.11) 




AJjJl 


ATj’ , , At? , . + Atj . . Atj . , 
i.J.k i,j,k i,j,k i,j,k 


+ At? 



for k « £, h, and t,. Because the eigenvalues are: computed at cell faces 

l 

rather than cell centers, is the average of the eigenvalues on cell 
faces in the k » constant computational planes (where k = n, or c) , 
and L is eigenvalue 4 or 5 because one of these will always have the 
maximum absolute value. To accelerate convergence for steady-state 
solutions the maximum allowable time step in each volume is used where 


CFL < 2. 


~n 



For the computation of the first points inside the computational 
domain, the scheme is only first order accurate and stable for a CFL 
of 1. The time steps at these outside points are, therefore, decreased 
by a factor of 2. 



viryr'V 'Ut •V***'* 
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IV. BOUNDARY CONDITIONS 


To derive the characteristic variable boundary conditions, the Euler 
equations need to be cast in characteristic variable form. Therefore, 
the characteristic variable form of the equations is derived first, 
followed by the derivation of the boundary conditions for the specific 
cases of supersonic, and subsonic inflow and outflow, and impermeable 
surfaces . 

4.1 Characteristic Variables 

Consider the nonconservative form of the Euler equations given 
by Eq. (2.6). Using Eq. (2.18) in Eq. (2.6), and multiplying Eq. (2.6) 
on the left by P ^ , gives 


p~^ la .i. p a p”i la + p“^ p j[ p”^ la s o 

*k Dt f 1 k k A k 1 k 3k *k m m r m 3m 


(4.1) 


where k *= £ , r), or Ct and m includes the remaining two curvilinear 
coordinates out of the set £, q, or £ where m 4 k. Therefore, there 
is a two-tenn summation on m in the third term in Eq. (4.1). Define 
the third term in Eq. (4.1) as 


S. * P' 1 ? A P" 1 
k,m lc mm m 3m 


(4.2) 


with the two-term summation on m understood. Equation (4.1) can be 
written 

p- 1 Is. 4 . A p"; 1 la + s =--• 0 (4.3) 

k k k 3k k,m 


-1 -1 

Consider P, to be a constant matrix denoted as P, . Then 
k k, o 

Eq. (4.3) can be written 
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-!- A, + S. “ 0 

3t k 3k k,ra 


(4.4) 


Defining the characteristic vector, W^, as 


W, = P7 1 q 
k k,o 1 


(4.5) 


Eq. (4.4) becomes 


3W 9W 

™ + A. -~ + S. =0 
3t k 9k k,m 


(4.6) 


The elements of the characteristic vector, W, , are called characteristic 

* k’ 


variables, and are denoted by w , . 

The charac 
The vector q is 


The characteristic variables, w, ., are determined using Eq. (4. 5), 

K f 1 


q “ dtp, m, v, w, p]' 


.-1 


and the matrix q is given by Eqs. (2.23) where the variables p and 
c (speed of sound) in Eq. (2.23) are denoted p and c c to indicate 
a reference condition. Using q and Pj^ 0 in Eq* (4.5), the elements 
of the characteristic vector 


W k (w k,l’ w k,2’ W k,3’ W k,4’ W k,5 3 


T 


(4.7) 


are 


w, B J — [k (p - P) + k v - k w-'J 

k , 1 i vk i - c 2 * y 


(4.8a ) 


w = [k (p - - k u + k w) 

’ |Vk| y c 2 o 


(4. 8b) 


w k 3 = T“7 f k 2 ( P - " P 2 ) + lc v U " k x v3 

k ’ j Vk | C o 7 


(4.8c). 
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|Vk| >'2 p o c o 


(k x u 


+ k 


+ k w) ] 


— i! i— [ _ (t u + It v + k »)) 

jVk( G p o C o X >' 2 


(A . 8d) 


(A . 8e) 


The characteristic variables correspond in order to the eigenvalues 


X, «ktt + kv + k»»ff, 
k x y st 1 

X k = V 


X k “ 6 k 
k k 


\ - \ + c!Vk| 
“ e k - c[vk| 


(A. 9a) 
(A. 9b) 
(A. 9c) 
(A. 9dt) 
(A. 9e) 


A. 2 Characteristic Variable Boundary Conditions 

The boundary conditions are derived below assuming locally one- 
dimensional flow. This assumption is probably better for far field 
boundary conditions than for boundary conditions applied on or near 
surfaces. However , numerical experiments using aero pressure gradient, 
extrapolation, normal pressure gradient, and locally one-dimensional 
characteristic variable boundary conditions,, indicate that similar results 
can be obtained vising any of these four methods for impermeable, wall 
boundary conditions as long as the grid is not extraordinarily coarse. 

For the computations performed thus far, the locally one-dimensional 
characteristic variable impermeable wall boundary conditions are to be 
preferred over the other three methods; hence, this method is developed 
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below. (Further investigations should be performed without the 
locally one -dimensional assumption, however.) 

By neglecting tire Jft directions in Eq. (4.6) one obtains 


aw, 

k 

9x 


+ A, 


aw, 


k 

9k 


0 


This equation can be written 


dW. 

1c 

dx 


aw, 

k 

3x 


+ A 


k 9k 


0 


(4,10) 


(4.11) 


where 



The eigenvalues, therefore, indicate a direction in computational 

space. According to Eq. (4.11) one particular eigenvalue is associated 

i 

with one. particular characteristic variable. Each eigenvalue, Aj_, 
indicates the direction across the k « constant computational surface 
that information contained in the associated characteristic variable, 
w, , propagates. This result is the basis for determining the 
boundary conditions referred to here as characteristic variable boundary 
conditions. 

Boundary conditions are now developed for the following five 
specific cases 

1. supersonic inflow 

2. supersonic outflow 

3. subsonic inflow 

4. subsonic outflow 


5. 


impermeable surface 



Supersonic Inflow 


This is a situation where all eigenvalues have the same sign. 

Because flow is coming into the computational domain, all flow variables 
are specified. 

Su personic Outflow 

This is another situation where all eigenvalues have the same sign. 
Because flow is leaving the computational domain all flow variables at 
the boundary must be obtained from the solution in the computational 
domain. All flow variables are extrapolated from inside the computational 
domain to the boundary. 

Sub s onic Inflow 

This situation is characterized by four eigenvalues of the same sign 
and one of differing sign. For the subsonic inflow case shown in Fig. 

3a with the flow in ilc direction of increasing computational coordinate 
k, the first four eigenvalues are positive and the fifth is negative. 

For the subsonic inflow case shown in Fig. 3b with the flow in the 
direction of decreasing computational coordinate k, the first three and 
fifth eigenvalues are negative and the fourth eigenvalue is positive. 

For a totally general three-dimensional code either situation in Fig. 3 
could occur. Each possiblity in Fig. 3 is talcen into account in the 
following derivation. Using the characteristic variables given by 
Eqs. (4.8) one obtains 


[k x (p 

- - 77 ) + k v.-'k w] *=* 
2 z y a 

c 7 

[k.(p - 

~2> + \V ~ k y W] b 
c J 

(4.12a) 


o 


O 


[k (p 

y 

- -|) - k z u + k x wj a - 
c 

[k y (p - 

-f) - t V .» + k.w] b 

C 

(4.12b) 


o o 



n<- 2 Cp - •*%>) + k u - k^v 3 a » [k z (p - - E ^> + k u - k^v]^ (4.12c) 


[ — I A (k u + k v + ic w)i = [—1 ^ ± (k u + k v + k w) }, (4.12d) 

l p c x y 2 ,J a l p c x y z b 

o o ' o o 7 


(k u •!• k v + k w))» = [-2^^- + (k u fk v + k w)L (4.12c) 
l p x y z ' J Z P„ c „ x y z b 


o o 


o o 


where the plus signs in Eq, (4. 12d) and the negative signs in Eqs. (4.12e) 
refer to the situation in Fig. 3a, and the other sign option refers to 
the situation in Fig. 3b. The subscript a refers to approaching the 
boundary, subsc.ipt b refers to the boundary, and subscript t refers to 
leaving the boundary (see Fig. 3). Note from the equations for k^, k^, 
and k given below Eq. (2.2) that the products Jk , Jk , and Jk„ are 
components of area vectors. The compute?: code actually uses these 
components of cell surface areas, and because the boundary of interest 
in Eqs. (4.12) is the cell face containing the boundary point b, the 
area vector components of interest correspond to this cell face. The 
metrics at points a and <£ are taken to be the same as those at. point b, 

«tnd the coefficients of the bracket terms in Eqs. (4.8) are not carried 
through the manipulation of the equations. The linearization point o 
is taken to bo on the boundary. 

Equations (4.12d) and (4..l2e) can he combined to obtain 


P b " ’“ {P a + * P o C o lk x (u a “ u £ } + k y (v a " v £ } + V w a “ U t ) ] } (4,13a) 


where k , k , and k are defined by Eqs. (2.15). Equations (4.12a), 
x y z 

(4.12b), (4,12c), and (4,12d) can be solved for the four remaining 
unknown boundary values, giving 
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P a “ P b 


u + k 


v ± k 


w ± k 
a ape 
o o 


p o C o 


P a ~ 

£b 

O 0 


P a " 

p b 


C$,1 3b) 


(4.13c) 


(4,13d) 


(4. 13e) 


The plus sign option in Eqs. (4.13) refers to the computational coor- 
dinate and inflow situation depicted in Fig. 3a, and the negative sign 
option refers to Fig. 3b. There is no sign option in Fig, (4,13b), 

Note that these signs correspond to the sign of the first three eigenvalues, 
and hence this is a means of writing the code for general applications 
with arbitrary orientation of the computation coordinates. The point 
a is outside the computational domain, point, b is on the computational 
boundary, and point t is inside the computational domain,. 

Subsonic Outflow 

This situation is also characterized by four eigenvalues of the 
same sign and one of opposite sign. The development of subsonic out- 
flow boundary conditions is similar to that for subsonic inflow, and 
Fig, 3 can be used again for illustration. However, for subsonic out- 
flow only one characteristic variable is specified and four are determined 
from information inside, the computational domain, whereas, for subsonic in- 
flow four characteristic variables were specified and one was determined from 
information inside the computational domain. Using the characteristic 
variables given by Eqs. (4. 8), and the signs of the eigenvalues given 
by Eqs. (4.9) for the situations in Fig. 3, one obtains a set of 
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equations identical to Eqs. (A. 12). Note, however, that although 
the equations are identical, point a is now inside the computational 
domain and point l is ouftside the computational domain; whereas, just 
the opposite was true for subsonic inflow. 

Because the equations for subsonic outflow are the same as Eqs. 


(A. 12),. they have the same formal solution. However, because one 
characteristic variable is specified, the resulting boundary conditions 
differ somewhat from the subsonic inflow boundary conditions. Consider 

Eq. (A.12e). By specifying that the outflow is straight, 
according to Eqs. (A.12e). The remaining four equations 
for the remaining four variables giving 

then P b “ P£ 
can be solved 

p b - p.e 

(A .lAa) 

, p b " Pa 

P b = p a + 2 

c 

o 

(A.IAb) 

if Pa ~ p k 
u b “ u a ik x pc' 
o o 

(A.IAc) 

~ P a “ p b 

v, « v ± k 

b a y p c 

(A.lAd) 

~ P a " P b 

w, = w ± k 

b a. z p c 

(A.lA e ) 


The plus and minus signs have the same meaning here as for Eqs. (A. 13). 
For external flow computations, p ^ could be the ambient static pressure, 

P . 
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Impermeable Surface 


\ 


For a boundary across which there is no flow the first three 


eigenvalues given by Eq, :^4>-€0 are zero, the fourth is positive, and 
the fifth is negative. One condition must, therefore, be specified. 
The condition specified is that there is no flow across the boundary. 
The following relations among the characteristic variables are used to 
determine the boundary conditions 


[k x (p 

1) + k v - k w] <* [k (p 

Z Z VO X 

c } 

- -|) + 
c 

fe v - k w] 

Z y r 

(4.15a) 


0 

o 



[k y (p 

~ *2 } - k z U + k x w] b = lk y (p 

.. -£-) - 
2' 

k u + k w] 
s x r 

(4.15b) 


c ' c 


[k (p 
*» 


-t) 

c o 


k u - k v], 
y X b 


o 

[k (P — ^r) +. & u - k v] 

2 2 y x r 

c 1 

o 


(A. 15 c) 


[k u + k v + k w], “ 0 
x y z b 


(4.15d) 


[-*>1 ? Ck „ 


P c 
o o 


k v 

y 


k z w)] b 


P I Vk 1 

P c ‘ 
o o 


+ (k u + 


k v + k w) ] 
y z r 


(4.15e) 


The subscript r refers to a reference value, which is selected as the 
center of the first cell from the boundary. The- minus and plus signs 
in Eq. (A.ISe) correspond to the location of the point r. If the point 
r is in the positive k direction from the boundary then the minus 
sign is used in Eq. (4.15e), and if it is in the minus direction then 
the plus sign is used. 
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Finite volume codes only require the pressure at an impermeable 
boundary and consequently Eqs. (4.15a), (4.15b), and (4.15c) are not 
needed. However, to Sie handling of points near boundaries 

and aid in code vectorization, phantom points are used in the present 
version of the code. The use of phantom points requires information of 
variables other than pressure to ensure, for example, zero flow across 
an impermeable boundary. Such information can be obtained from Eqs. 
(4.13). 

Equations (4.15d) and (4.15e) can be solved for p^. Equations 
(4.15a) - ( 4.15d) can then be solved for the remaining four variables. 
The solution of Eqs. (4.15) is 


p, ■= p +pc(ku +kv -i-kw) 
b r o o x r y r z r 


P, - P 
b r 


P b - P r + 2 


(4.16a) 

(4.16b) 


u, = u - k (t u + k v + k w ) 
b r x'xr yr z r 


(4 .16c) 


v « v - k (k u +kv *tkw) 
b r y x r y r z r 


(4 . 16d) 


w, K w •- k (k u + k v + k w ) 
b r z xr yr zr 


(4. 16e) 


where the point r is the center of the first cell from the boundary 
and the minus sign in Eq. (4.16a) is used if r is in the positive k 
direction from the boundary, and the plus sign is used if r is in the 
negative k direction from the boundary. 
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A. 3 Phantom Points 

Phantom points are denoted by the subscript p. The points are 
obtained from the relations 

P p “ 2P b " P in (4 * 17a) 

p p " 2p b " p in (4 ' 17b) 

U p = 2u b - U in (4 * 1?C) 

v P =2v b" v in (4 - 17d) 

w *» 2w u - w, (A.17e) 

p b in 

where the subscript in refers to the center of the first cell inside 
the computational domain and can be any of the points a, ■£, or r used 
in this section. For example, the phantom points for an impermeable 
surface are 


p = 


+ 2p c (k u 

+ k v 

k w ) 

(4.18a) 

p 

r 

o o x r 

y r 

z r 




2(P b “ P r ) 




p„ B 


+ S—JL. 



(4.10b) 

p 

r 

2 

c 






o 




u ■= 

u 

- 2k (k u + 

k v + 

k. w ) 

(4.18c) 

p 

r 

x x r 

y r 

z r 


V a 

V 

- 2k (k u + 

k v + 

k w ) 

(4.18d) 

P 

• r 

y x r 

y r 

z r 


W a 

w 

- 2k (k u + 

k v + 

kw) 

(4.18e) 

p 

r 

z x r 

y * 

z r 



The velocity vector components are the same as those used by Jacocks and 
Kneile (Ref. 14). 
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V. RESULTS 

The computer program written to solve the three-dimensional un- 
steady Euler equations is referred to as the SKOAL code. SKOAL is 
an acronym for nothing. 

To investigate the results of this code a numerical solution was 
obtained for the ONERA M6 wing at = O .84 and a « 3.06°. The pressure 
distribution is shown in Fig. 4, This solution is compared to a 
solution from the FL057 code in Fig. 5. Identical 96x16x16 meshes 
were used for both solutions in Fig. 5. The major difference between 
the two solutions occurs in the outboard region of the wing. More 
detailed comparisons between the solutions , including comparisons 
with experimental data (Ref. 15), are given in Fig. 6. Comparisons of 
spanwise distributions of lift and drag are given in Fig. 7. The 
FL057 code gives a 4% higher value of lift to drag ratio than the 
SKOAL code. 

The number of supersonic points (NSUP) Is sometimes used as a 
crude indication of convergence. The FL057 solution was obtained by 
first using a 48x8x8 grid and then using this crude grid solution as 
initial conditions for the 96x16x16 grid solution;. Hence the NSUP 
corresponding to an impulsive start for the 96x16x16 grid was not 
available from the FL057 code for comparison with the SKOAL code. In 
order to compare the NSUP history between the two codes, a 48x8x8 grid 
solution was obtained using the SKOAL code. These results are presented 
in Fig. 8. The FL057 code is stable for a CFL of 2.0 using a four- 
stage Runge-Kutta scheme, whereas, the SKOAL code is stable for a CFL 
of 2.0 using the predictor-corrector scheme. It is interesting to note 
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that the number of cycles required for the NSUP to become constant 
using FL057 is 246, whereas, the number of cycles required for the 
NSUP to become constant using SKOAL is 325. The ratio of these two 
numbers is essentially the Inverse ratio of CFL numbers. Another way 
of looking at this is to note that FL057 passes through a flux balance 
routine four times during each cycle, whereas SKOAL passes through a 
flux balance routine three time during each cycle. A comparison of 
the ratio of the NSUP to the final number of NSUP (denoted as NSUP/ 
(NSUP) c ) as a function of flux balances is given in Fig. 9» Based on 
the number of flux balances required to reach steady state for this 
solution these methods are operating identically. However, one pass 
through a dissipation routine (which is similar in terms of computer 
resources to an extra pass through the flux balance routine) is required 
in the FL057 code wnich is not required in the SKOAL code because no 
additional dissipation terms were included. The important comparison, 
however, is the number of computer resource units required to reach 
steady state. This depends on the coding, vectorization, storage, etc., 
of each code on the same machine. Such a comparison cannot presently 
be made. As the codes now stand, the FL057 code is probably at least 
twice as fast as the SKOAL code per cycle. However, based on Fig. 9, 
it is anticipated that the SKOAL code can be improved. 

Numerous other results have been obtained using the SKOAL code 
for various two-and three-dimensional geometries, and for subsonic, 
transonic, and supersonic flow. The ONERA wing, however, is the only 
solution obtained thus far that can be compared to another Euler code 
solution for which the same, mesh was used. 



VI. CONCLUDING REMARKS 


A method was presented for solving the three-dimensional unsteady 
Euler equations, bas,ed on flux vector splitting. The equations were 
cast in curvilinear coordinates and a finite volume discretization 
was used for handling arbitrary geometries. The discretized equations 
were solved using an explicit, upwind second-order predictor-corrector 
scheme that is stable for a CFL of 2. No additional dissipation terms 
were included in the scheme. Local time stepping was used to accelerate 
convergence forsteady-state problems. Characteristic variable boundary 
conditions were developed and applied in the far field and at impermeable 
surfaces. Numerical results were obtained and compared with results 
from the FL057 code and experimental data for the ONERA H6 wing at 
M =0.84 and cs = 

oo 
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Fig. 3 Computational Coordinates Schematic for Characteristic 
Variable Boundary Conditions. 
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Numerical and Experimental Pressure Distributions at Various Span Locations. 
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Fig. 7 Bpanwise Lift and Drag Distribution. 
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Fig. 8. Number of Supersonic Points (NSUP) as a Function of 
Cycles. 





QNERA MS WHIG 
ft * om 




Fig.. 9 Number ox Supersonic Points (NSUP) as a Function of Flux 
Balances . 
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